rm(list=ls())
gc()

library(dplyr)
library(rstan)
options(mc.cores = parallel::detectCores())
######### 1. Run IRT (Time lag = year or legislature) ########


make_IRT_data <- function() {
  
  swissvotes <- read.csv(url("https://swissvotes.ch/page/dataset/swissvotes_dataset.csv"), sep = ";")
  
  colnames(swissvotes)
  swissvotes <- swissvotes[swissvotes$anr>=193 & ### After 1965
                             swissvotes$anr<=664 & ###Before 2015 
                             swissvotes$rechtsform<=3, ## Only initivative, referenda mandatory and facultative
                           c(which(colnames(swissvotes)=="anr"),
                             which(colnames(swissvotes)=="datum"),
                             which(colnames(swissvotes)=="p.fdp"),
                             which(colnames(swissvotes)=="p.cvp"),
                             which(colnames(swissvotes)=="p.svp"),
                             which(colnames(swissvotes)=="p.gps"),
                             which(colnames(swissvotes)=="p.sps"),
                             which(colnames(swissvotes)=="rechtsform"), 
                             which(colnames(swissvotes)=="legislatur"))]
  
  for (i in 3:7) {
    if (i==3) {
      data <- swissvotes[,c(1, 2, 9)]
      data$pos <- swissvotes[,i]
      data$party <- substr(colnames(swissvotes)[i], 3, 5)
      data$year <- as.numeric(substr(swissvotes$datum, 7, 10))
    } else {
      dta <- swissvotes[,c(1, 2, 9)]
      dta$pos <- swissvotes[,i]
      dta$party <- substr(colnames(swissvotes)[i], 3, 5)
      dta$year <- as.numeric(substr(swissvotes$datum, 7, 10))
      data <- rbind(data, dta)
    }
  }
  
  
  data <- data %>% 
    mutate(i = as.numeric(as.factor(anr)),
           j = as.numeric(as.factor(party)),
           t = as.numeric(as.factor(year)),
           t2 = as.numeric(as.factor(legislatur)),
           y = ifelse(pos==1, 1, ifelse(pos==2, 0, NA)))
  
  y <- data$y
  NAS <- which(is.na(y))
  y <- y[-NAS]
  i <- data[-NAS,]$i
  j <- data[-NAS,]$j
  t <- data[-NAS,]$t
  t2 <- data[-NAS,]$t2
  N <- length(y)
  I <- max(i)
  J <- max(j)
  T <- max(t)
  T2 <- max(t2)
  
  data_IRT <- list(N=N, I=I, J=J, T=T, 
                   y=y, i=i, j=j, t=t)
  data_IRT2 <- list(N=N, I=I, J=J, T=T2, 
                    y=y, i=i, j=j, t=t2)
  return(list(data_IRT=data_IRT, data_IRT2 = data_IRT2))
}


data_IRT <- make_IRT_data()

### 2. Run the IRT model ####

run_IRT <- function(data_IRT) {
  
  model <- stan(file = "Dynamic IRT.stan", data = data_IRT, iter = 10000, warmup = 5000,
                seed = 12345, cores=4, thin = 1)
  
  return(model)
  
}


run_IRT_simple <- function(data_IRT) {
  
  model <- stan(file = "Simple IRT.stan", data = data_IRT, iter = 10000, warmup = 5000,
                seed = 12345, cores=4, thin = 1)
  
  return(model)
  
}


model1 <- run_IRT_simple(data_IRT$data_IRT2)

model2 <- run_IRT(data_IRT$data_IRT2)


save(model1, file = "model/Model party ideal position.rda")
save(model2, file = "model/Model party ideal position_by legislature.rda")






